function [result_beam, result_range, result_refFreq] = GetCellAbsorption(obj, beamName, varargin)
    p=inputParser;
    p.addParameter('frequencyRange', [], @(x) isempty(x) || isnumeric(x));
    p.addParameter('isPlot', false, @islogical);
    p.parse(varargin{:});
    
    
    try
        beam = obj.system.getCompomentByName(beamName).stuff;
    catch
        error('Cannot obtain beam object with name = [%s]', beamName);
    end
    
    result_beam = calc_absorption(obj, beam.parameter.detuning, beam.refTransition);
    result_refFreq = beam.refAtom.parameter.reference_frequency(beam.refTransition);
    
    if ~isempty(p.Results.frequencyRange)
        result_range = calc_absorption(obj, p.Results.frequencyRange, beam.refTransition);
    else
        result_range = [];
    end
    
    
    if p.Results.isPlot && ~isempty(p.Results.frequencyRange)
        f = figure;

        f0 = result_beam.frequency/1e3;
        cs_beam_average = result_beam.cross_section_average;
        od_beam_average = result_beam.optical_depth_average;
        ts_beam_average = result_beam.transmission_average;
        
                
        ax1 = subplot(3, 1, 1, 'Parent', f);
        plot(ax1, result_range.frequency/1e3, result_range.cross_section, '--'); grid(ax1, 'on'); hold(ax1, 'on');
        plot(ax1, result_range.frequency/1e3, result_range.cross_section_average, 'k-', 'linewidth', 1); 
        plot(ax1, [f0, f0], [0 1.2*max(result_range.cross_section_average)], 'k--', f0, cs_beam_average, 'kd', 'MarkerFaceColor', 'black');
        hold(ax1, 'off');
        ylabel(ax1, 'cross section (cm^2)');
        legend(ax1, [result_range.names, 'average'], 'Location', 'best');
        title(ax1, sprintf('cross section = %3.2e cm^2', cs_beam_average));
        
        ax2 = subplot(3, 1, 2, 'Parent', f);
        plot(ax2, result_range.frequency/1e3, result_range.optical_depth, '--'); grid(ax2, 'on'); hold(ax2, 'on');
        plot(ax2, result_range.frequency/1e3, result_range.optical_depth_average, 'k-', 'linewidth', 1); 
        plot(ax2, [f0, f0], [0 1.2*max(result_range.optical_depth_average)], 'k--', f0, od_beam_average, 'kd', 'MarkerFaceColor', 'black');
        ylabel(ax2, 'optical depth');
        hold(ax2, 'off'); legend(ax2, [result_range.names, 'total'], 'Location', 'best');
        title(ax2, sprintf('OD = %3.2f', od_beam_average));
        
        ax3 = subplot(3, 1, 3, 'Parent', f);
        plot(ax3, result_range.frequency/1e3, result_range.transmission, '--'); grid(ax3, 'on'); hold(ax3, 'on');
        plot(ax3, result_range.frequency/1e3, result_range.transmission_average, 'k-', 'linewidth', 1); 
        plot(ax3, [f0, f0], [0 1], 'k--', f0, ts_beam_average, 'kd', 'MarkerFaceColor', 'black');
        hold(ax3, 'off');
        xlabel(ax3, sprintf('Detuning (GHz) [%9.6f THz]', result_refFreq/1e12)); ylabel(ax3, 'transmission');
        legend(ax3, [result_range.names, 'total'], 'Location', 'best');
        title(ax3, sprintf('Transmission = %3.2e', ts_beam_average));
    end
end

function result = calc_absorption(obj, frequency, transition)
    abundanceList = arrayfun(@(v) v.abundance,                            obj.system.getVaporList());
    densityList   = arrayfun(@(v) v.density*1e-6,                         obj.system.getVaporList()); % cm-3
    result.names  = arrayfun(@(v) v.name,                                 obj.system.getVaporList(), 'UniformOutput', false);    
    crs_section   = arrayfun(@(v) v.cross_section(frequency, transition), obj.system.getVaporList(), 'UniformOutput', false);
    
    result.frequency = frequency;
    result.cross_section = cell2mat(crs_section);
    result.cross_section_average = result.cross_section*abundanceList';
    
    result.optical_depth = result.cross_section .* densityList * obj.conditions.cellSize;
    result.optical_depth_average = sum(result.optical_depth, 2);
    
    result.transmission = exp(-result.optical_depth);
    result.transmission_average = exp(-result.optical_depth_average);
end